% Written to plot resitivity/phase patches (sections)
% To be used after calling Prof = Resprof.m 
% latest adet 1st sep 2004

mode = 2; % 1 res 2 phase

colormap('default'); % because prior changes may affect the color codings
%C = colormap;
nstn = length(Prof);
clear Data;
if mode == 1,
    for i = 1:nstn,
    Data(i,:) = Prof(i).arr_res;
    end;
    Data = log10(Data);
    mindata = log10(10);
    maxdata = log10(100000);
    disp('Enter A title string');
    tit = input([Prof(1).DataCode '_Res :']);
  elseif mode == 2,
    for i = 1:nstn,
    Data(i,:) = Prof(i).arr_ph;
    end;
    mindata = 0;
    maxdata = 90;
    disp('Enter A title string');
    tit = input([Prof(1).DataCode '_Phase :']);
end;


for i = 1:nstn,
X(i) = Prof(i).dist;
end;

%Y = log10(1./Prof(1).master); % Y is log10(period) in second
Y = log10(Prof(1).master);

xx = [X(1)-(X(2)-X(1))/2 X(1:end-1)+diff(X)/2 X(end)+(X(end)-X(end-1))/2];
yy = [Y(1)-(Y(2)-Y(1))/2 Y(1:end-1)+diff(Y)/2 Y(end)+(Y(end)-Y(end-1))/2];

C = flipud(colormap); % Reverses -> red for conductor and blue for resistor
kk=1;
for i = 1:length(X),
for j = 1:length(Y),
if isnan(Data(i,j)) ~= 1, 
    Xpatch(kk,:) = [xx(i) xx(i+1) xx(i+1) xx(i)];
    Ypatch(kk,:) = [yy(j) yy(j) yy(j+1) yy(j+1)];
    cidx = floor(63*(Data(i,j)-mindata)/(maxdata-mindata)+1);
    
    if cidx > 64,
        cidx = 64;
    elseif cidx < 1,
        cidx = 1;
    end;
    Cpatch(1,kk,:) = C(cidx,:);
    kk=kk+1;
end;
end;
end;


for i = 1:length(X),
%SS1(i,:) = [Prof(i).file(1:2) Prof(i).file(5:6)]; % this may have to be changed for file names with different nomencl.
SS1(i,:) = [Prof(i).file(1:4)];
end;

h=patch(Xpatch',Ypatch',Cpatch);
set(h,'Edgecolor','none');
axis([-5.0000  180.0000   -2.2    3])
a = axis;
h=text(X,(a(end)-0.6)*ones([1,length(X)]),SS1);
set(h,'Rotation',90);
colormap(C);
caxis([mindata maxdata]);
colorbar('horiz');

ylabel('Log10(Frequency) Hz');
xlabel('Distance (km)');
title(tit);

% hold on;
% for i = 1:length(X),
%     plot(ones([1,length(Prof(i).frq)])*X(i),log10(Prof(i).frq),'r*');
% end;    
